Loading the coordinates data

We need the sf package:

> if (! "sf" %in% rownames(installed.packages())) install.packages("sf")
> library(sf)
Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3

The coordinates of the houses of the Hanam cohort are in a zipping shapefile that can be downloaded as so:

> download.file("https://www.dropbox.com/s/yyg6gx0nycpv60b/HH%20point%20Hanam.zip?raw=1", "hanam.zip")

Once downloaded, unzip the file:

> unzip("hanam.zip")

Once unzipped read the shapefile:

> hanam <- sf::st_read("Household_point.shp")[, "Name"]
Reading layer `Household_point' from data source `/Users/choisy/Dropbox/aaa/projects/Hanam/Household_point.shp' using driver `ESRI Shapefile'
Simple feature collection with 263 features and 22 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 105.9174 ymin: 20.48691 xmax: 105.9356 ymax: 20.51665
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

And clean the disk

> file.remove(grep("^Household", dir(), value = TRUE))
[1] TRUE TRUE TRUE TRUE TRUE TRUE

and

> file.remove("hanam.zip")
[1] TRUE

Here the data, which is of sf class

> hanam
Simple feature collection with 263 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 105.9174 ymin: 20.48691 xmax: 105.9356 ymax: 20.51665
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
   Name                  geometry
1  H001 POINT (105.9274 20.49189)
2  H002 POINT (105.9286 20.49196)
3  H003 POINT (105.9314 20.49663)
4  H004 POINT (105.9274 20.49191)
5  H005 POINT (105.9303 20.49529)
6  H006  POINT (105.9329 20.4953)
7  H007 POINT (105.9296 20.49722)
8  H009 POINT (105.9266 20.49759)
9  H010 POINT (105.9296 20.49672)
10 H011 POINT (105.9306 20.49684)

The coordinates of the houses can be extracted from that object and converted into a data frame as so:

> houses <- cbind(house = hanam$Name, as.data.frame(st_coordinates(hanam)))

which gives:

> head(houses)
  house        X        Y
1  H001 105.9274 20.49189
2  H002 105.9286 20.49196
3  H003 105.9314 20.49663
4  H004 105.9274 20.49191
5  H005 105.9303 20.49529
6  H006 105.9329 20.49530

Plotting the houses

We need the OpenStreetMap package:

> if (! "sf" %in% rownames(installed.packages())) install.packages("OpenStreetMap")
> library(OpenStreetMap)
Registered S3 methods overwritten by 'ggplot2':
  method         from 
  [.quosures     rlang
  c.quosures     rlang
  print.quosures rlang

Below is the code you’d need to download the tiles from 3 different types of maps:

> upperleft  <- c(20.525681, 105.905383)
> lowerright <- c(20.463343, 105.969019)
> bing <- openmap(upperleft, lowerright, type = "bing", minNumTiles = 20)
> osm  <- openmap(upperleft, lowerright, type = "osm",  minNumTiles = 20)
> esri <- openmap(upperleft, lowerright, type = "esri", minNumTiles = 20)

But the downloading would be twice as fast if you downlaod the objects directly from here:

> download.file("https://www.dropbox.com/s/26y2pgouodj4rpe/bing.rds?raw=1", "bing_file.rds")
> download.file("https://www.dropbox.com/s/rd34k3ixwthzg9g/esri.rds?raw=1", "esri_file.rds")
> download.file("https://www.dropbox.com/s/827ze3xr0orbab0/osm.rds?raw=1" ,  "osm_file.rds")
> bing <- readRDS("BING.rds")
> esri <- readRDS("ESRI.rds")
> osm  <- readRDS("OSM.rds")
> for(file in paste0(c("bing", "esri", "osm"), "_file.rds")) file.remove(file)

The function that makes the map:

> plot_hh <- function(map, points) {
+   require(sf)
+   plot(st_geometry(st_transform(points, map$tiles[[1]]$projection)))
+   plot(map, add = TRUE)
+   plot(st_geometry(st_transform(points, map$tiles[[1]]$projection)), add = TRUE, col = "red")
+ }

OSM map:

> plot_hh(osm, hanam)

ESRI map:

> plot_hh(esri, hanam)

BING map:

> plot_hh(bing, hanam)

Calculating distances

You can calcultate distances between 2 points with the sf::st_distance(). For that, you need to install the lwgeom package:

> if (! "lwgeom" %in% rownames(installed.packages())) install.packages("lwgeom")

Once this is done, the calculation of the distance between the first 2 houses is done as so:

> sf::st_distance(hanam[1, ], hanam[2, ])
Units: [m]
         [,1]
[1,] 129.5865

You can get the distance matrix directly like so:

> distmat <- sf::st_distance(hanam)

Clustering:

> hc <- hclust(as.dist(distmat))

Then, you can generate clusters, depending on a distance cut-off. For example with a cut-off of 1000 m, your clusters would be:

> cutree(hc, h = 1000)
  [1] 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 2 1 1 1 2 1 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 2 1 2 1 1 1 1
 [67] 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 4 4 4 4 4 4 4 4 4 4 4 4 2 1 4 4 4 4 4 4 4 5 3 3 6 6 3 3 3 3 3 3 6 3
[133] 6 6 3 5 3 6 6 3 6 3 6 6 5 5 6 6 6 3 6 3 6 3 3 6 6 3 6 5 3 3 3 3 3 6 3 3 3 3 3 6 3 3 6 6 3 3 3 6 6 6 6 6 6 6 6 6 6 7 7 8 8 8 8 8 8 8
[199] 3 8 8 7 8 7 8 5 8 8 8 3 8 8 8 8 8 8 7 7 7 7 7 7 5 8 7 8 9 9 9 1 1 1 3 9 6 6 6 9 6 3 1 5 3 6 3 9 3 6 9 9 9 9 9 5 1 5 1 2 1 1 6 6 6

You can merge this information with the house identifier in a data frame as so:

> clusters <- data.frame(house = hanam$Name, cluster = cutree(hc, h = 1000))

Which gives:

> head(clusters)
  house cluster
1  H001       1
2  H002       1
3  H003       1
4  H004       1
5  H005       1
6  H006       2